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[DESCRIPTION] z 

FIELD OF THE INVENTION 

5 

The present invention relates to a method for detecting the 
orientation of object shapes in radiographic images. 

BACKGROUND OF THE INVENTION 

10 

In radiological practice, it is common to display different 
exposures pertaining to a patient study in a predefined format. This 
feature is known as a hanging protocol. In film-based operation, it 
means that the radiologist or operator hangs the films on the light 
is box in a specific spatial arrangement according to standard or local 
preferences. Determination of the orientation or reflection of an 
examination type, or verification of it when it is available, is 
beneficial to the correct display of many examination types. 



20 For example, thorax examinations are always hung in such a way 

as to display the ribcage in the upstanding direction and showing 
the heart shadow at the right side with respect to the centerline. 
However, for reasons of exposure and read-out geometry, 90 -degree 
orientation changes may occur in digitally acquired images. 

25 In particular for chest images acquired with a computed 

radiography modality, guaranteeing a correct orientation is a major 
problem, because the cassette may be wrongly oriented at exposure 
time due to patient conditions and clinical requirements. 
Orientations different from the standard upright orientation are 

30 reported to be as high as 40%. Especially in high throughput 

hospital departments, radiologist and clinicians perceive the task 
of rotating the image to the standard orientation as a cumbersome 
task and a waste of time. Moreover, when automatic diagnostic data 
processing has to be performed on the digital image, a standard 

35 orientation must be available for correct operation of the CAD 

(computer-aided diagnosis) algorithms or subsequent operations such 
as annotations, measurements and archiving. 
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A similar problem exists in screening mammography, where 
standard 4 -views of left and right breast are taken for each of the 
incidence directions (cranio- caudal (CC) and medio- lateral oblique 

5 (MLO) ) . These views are typically displayed in a mirrored fashion, 
such that the thorax edges or both breasts are central and touching, 
and the left breast image being displayed on the right and the right 
breast image being displayed on the left. However, because both 
breasts images are acquired in a similar manner, and it is in 

w general not known which image is either corresponding to the left or 
the right breast, one image must be flipped before it can be 
positioned adjacent to the other image. In conventional screen-film 
imaging, X-ray opaque lead letters are radiographed simultaneously 
(RCC, LCC, RMLO and LMLO) , and the RCC resp. RMLO films are flipped 

15 manually prior to hanging them on the right of the LCC resp. LMLO on 
the light box. 

Digitally acquired mammograms may still be read in a 
conventional way by printing them on film and displaying them on a 

20 light box. Pairs of mammograms (e.g. the RCC/LCC pair and the 

RMLO/LMLO pair) may be printed on a single large film sheet or on 
two smaller sized sheets. Generally, the print margin of a hardcopy 
machine is adjustable, so as to minimize the non-printed portion of 
an image. For mammography hardcopy, the print margin corresponding 

25 to the thorax side is kept as small as possible, so that a right- 
left pair of images, when viewed simultaneous and in close vicinity, 
shows a minimal non-diagnostic viewing area in between both images. 
Therefore, when using a pair of small sheets to print left and right 
image respectively, means to identify the thorax side automatically 

30 prior to printing them is needed, because the thorax side position 
is generally not known or it may not be assumed to be known. 
Likewise, in the large film option, where both images are printed on 
one file sheet in order to compose the image such that the right 
image is touching in a mirroring manner to the left image, knowledge 

35 of the thorax side of both images is needed as well. 
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Digital mammography may be read on a computer display or 
viewing station without having resort to printed mammograms, a 
viewing condition known as softcopy reading. However, also here, the 
sub- examination types identifying right and left images may not be 
5 known at display time. Furthermore, the thorax edge orientation may 
not be standardized, e.g. it may either touch the left or right or 
the upper or lower image border. Hence, there is a need to 
accomplish the mirrored viewing disposition in an automated way. 

io Generally, a hanging protocol function enables users to arrange 

and display images on medical viewing stations to suit specific 
viewing preferences. To this purpose, the sub -examination is used to 
assign the sub-images pertaining to a patient study of a body part 
to a position in a preferred display layout of that examination 

j 5 type. When the image sub- type is known, hence its position in the 

layout is determined, the image can still by oriented in 8 different 
ways: it can be oriented correctly, or it can be rotated of 90, 180 
or 270 degrees; in any of these four cases, the image can also be 
flipped around a vertical (or horizontal axis) . Therefore, there is 

20 a need to derive the orientation of the image automatically, to 
assure viewing according to the radiological standard or local 
preferences . 

The problem of orientation detection and occasional correction 
25 therefore consist in detecting the orientation as one of 0 , 90, 180 
or 270 degrees with respect to an image coordinate axis. When the 
orientation equals 0 degrees, no orientation correction needs be 
performed, otherwise the image must be anti-rotated to have it in 
the upright direction (head-neck region on top of the corrected 
30 image) . 

The orientation detection solution must further be insensitive 
to a number of disturbing image characteristics such as 
the presence of collimation borders 
35 - rotation of the centerline of the thorax with respect to the 

vertical axis 
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- the presence of other non-thorax regions, such as head/neck, 
sub-diaphragm region and abdomen region 

- unilateral hyper-dense appearance of a diseased lung region 

5 

SUMMARY OF THE INVENTION 

x The above-mentioned objects are realised by a method as set out 
in claim 1 . 

w 

Further embodiments are set out in the dependent claims. 

The method is particularly useful for orienting and 
standardizing chest radiographs in up-right position. 

15 

The method of the current invention is particularly useful in 
the field of digital mammography, to detect the thorax to nipple 
direction and align a left-right pair of breast views such that 
their thorax sides are touching when printed on film or displayed on 
20 screen. 

In the field of general radiology, the invention is 
particularly useful to determine the orientation of curved shaped 
objects, such as determination of the orientation of a frontal or 
25 lateral exposure of the skull in 2 D or 3D cerebral examination (MR, 
CT) . 

The embodiments of the methods of the present invention are 
generally implemented in the form of a computer program product 
so adapted to carry out the method steps of the present invention when 
run on a computer. The computer program product is commonly stored 
in a computer readable carrier medium such as a CD-ROM. 

Alternatively the computer program product takes the form of an 
electric signal and can be communicated to a user through electronic 
35 communication. 

BRIEF DESCRIPTION OF THE DRAWINGS 
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Fig. 1 shows a general overview of the method according to the 
present invention, 

Fig. 2 illustrates the orientation computed by edge gradient, 
5 Fig. 3 shows a planar curve, with local first and second derivative. 
Two coordinate systems are constructed at point pq : (1) the Frenet- 
Serret frame, having axes in the direction of the tangent to the 
curve, and perpendicular to the tangent direction (towards the 
center of the osculating circle) . (2) a local analyzing coordinate 
w frame, which is aligned with the Cartesian image coordinate axes in 
this example . 

Fig. 4 (a) and (b) illustrate the calculation of the curvature (a) by 
means of local subtended angle and (b) by means of the magnitude of 
the addition vector of 2 chord vectors, 
15 Fig. 5 is an example of orientation computation for thorax 
examination, 

Fig. 6 (a) to (d) illustrate an example of orientation computation 
for a mammographic image . 

20 DETAILED DESCRIPTION 

A general overview of the system to manipulate a radiation 
image according to the shape of objects in the diagnostic region (s) 
of the image is given in fig. 1. 

25 

In order to conduct a topological analysis of a radiographic 
image to infer the position of the breast thorax side and the 
nipple- thorax direction, shape analysis is performed. 

30 Shape analysis can be carried out starting from an intermediate 

representation typically involving the segmented image and/or 
special shape descriptors. 

A radiation image typically consist of three different areas: 
35 - The diagnostic area comprises pixels corresponding to patient 

anatomy. In general, the outline of this imaged area may take any 
shape. 
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- The direct exposure area is the image region that has received 
unattenuated radiation. Although this region has constant 
intensity corrupted by noise only, inhomogenit ies in incident 
energy (e.g. X-ray source Heel effect) and receptor (e.g. varying 

5 storage phosphor sensitivity in computed radiography) may distort 

this pattern. In European patent application EP 1 256 907 a method 
is disclosed to estimate these global inhomogenit ies 
retrospectively from the diagnostic image and flatten the response 
in all image parts in accordance with an extrapolated background 

10 signal. 

- The collimated areas appear on the image as highly attenuated 
pixels. The shape of these areas typically is rectilinear, but 
circular or curved collimation shapes may be applied as well. 

is Three different area transition types may be considered in a 

radiation image: diagnostic/direct exposure, diagnostic/collimated 
area, and direct exposure /coll imat ed area boundaries. 

Segmentation algorithms aim at detecting and separating of the 

20 set of pixels that constitute the object (s) under analysis. These 
techniques may be broadly classified according to the type of 
processing applied to the image. Region-based algorithms group 
pixels in the image according to suitable similarity criteria. In 
European patent application EP 887 769 a region-based algorithm is 

25 disclosed to segment direct exposure areas by grouping pixels 

according to centroid clustering of the gray value histogram. Edge 
based algorithms separate image pixels in high contrast regions in 
the image according to gray value differences of neighboring 
regions. In European patent application 610 605 and European patent 

30 application 742 536 an edge-based algorithm is disclosed to detect 

and delineate the boundaries between collimated areas and diagnostic 
areas on a single or multiply exposed image. Either in region-based 
and edge-based approaches, models may be used to restrict the 
appearance or shape of the segmented image areas to obey predefined 

35 photometric or geometric constraints. An example of this paradigm 

are the so-called Active Appearance and Active Shape Models (AAM and 
ASM) . 
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As a consequence of the segmentation result being either a 
region or a region transition, shape analysis techniques may also be 
divided roughly in either region-based or edge based procedures. 
5 Shape analysis techniques generally depart from a suitable 

representation of regions and region boundary, and hence they may 
also broadly be divided in region-based and contour-based classes. 
Examples of representations of either type are given in the sequel . 

io Shape analysis techniques are typically selected in view of the 

problem requirements. These requirements will broadly fall in one of 
two classes. The first class may be termed topology problems, in 
that the problem is one of locating and substantiating the specific 
spatial arrangement of an object shape with respect to other objects 

15 of reference systems. The second class may be termed 

characterization, and is closely linked with classification. In this 
problem class, shape specific descriptors must be applied, and in 
general a topology- independent characterization is needed. Both 
problem classes are closely linked, because when the topology of a 

20 specific shape needs be calculated, the shape must first be searched 
for and detected on the basis of the shape's specific 
characteristics. Conversely, when the characteristics of a specific 
shape need be determined e.g. in order to classify it, the shape 
must first be located in the image, which problem component is one 

25 of topology. Given the topology and characteristics of the shape, 
the application specific problem can be solved. 

The result of a shape analysis is a set of shape descriptors 
that characterize either the topology, the specific shape features 
30 or both. 

For the application of aligning a shape to a standard topology, 
the geometric transformation is determined by the geometric 
parameters of the shape's topology and the desired topology, and 
35 techniques known in the prior art for determining the geometric 
displacement field and the intensity interpolation are used to 
geometrically modify the image to the target topology. 



r 
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In order to determine the orientation of an object in a medical 
image, e.g. the thorax side of the breast mass, shape analysis 
techniques are used to describe the topology and characteristics of 

5 the object in the image. 

The shape of an object determines the extent of the object in 
the image, which is a binary image, and the spatial distribution of 
gray levels inside the object's extent. A shape analysis therefore 
departs from a shape representation, from which shape descriptors 

io are calculated. Shape representation methods result in a non-numeric 
representation of the original shape capturing the important 
characteristics of the shape (importance depending on the 
application) . Shape description methods refer to methods that result 
in a numeric descriptor of the shape, generated by calculating a 

15 shape descriptor vector (also called a feature vector) . The goal of 
description is to uniquely characterize the shape using its 
descriptor vector, independent of the shape's position, orientation 
and size. Conversely, the process of reduction of a shape to its 
canonical form that is invariant to translation, rotation and scale, 

20 assumes that the actual shape's position, orientation and size are 
determined . 



25 

Curve -based shape orientation measures 

Curve representation 

A curve or contour may in its simplest form be represented by 
30 the set of (possibly chained) contour pixels. At a higher level, the 
curve may be approximated in primitive forms such as a collection of 
approximating line segments (a polygonal representation, 
alternatively represented by the corner intersections), circle arcs, 
elliptical arcs, syntactic primitives, B-splines, Snakes and active 
35 contours, or multiscale primitives. Finally, a curve may be 

represented parametrically , for example as a two-component vector 
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y{t) = {x(t\y{t)} for a plane curve, or as a complex signal 
u (t) = x(t) + jy(i) , a chain code or a run- length code. 

Edge gradient measures 

5 For a continuous image f(x,y) , its derivative assumes a local 

maximum in the direction of the edge. Therefore, one edge detection 
technique is to measure the gradient of / along a straight line 
making an angle 6 w.r.t. the x-axis. If r is a parameter denoting 
arc length along the line, then the change df when a change dr is made 

io along the line is (see fig. 2) 

V = Vdx + ¥dy = fxCOS0+fin0 

dr • dx dr dy dr 

The maximum value in gradient df I dr is obtained by solving 
9 from 



15 which maximum value is 



\drj max v 



0 g ( x >y) ^ s the edge or gradient orientation, associated with the 

edge or gradient magnitude g(x,y) at position (x, y). 

Based on these concepts of continuous edge gradient, two types 

20 of edge detection operators for digital images have been developed. 
The first are gradient operators, based on calculations of the 
partial derivatives in x and y direction using finite differences. 
The Roberts, Prewitt, Sobel and isotropic operators compute 
horizontal and vertical differences of local sums to reduce the 

25 effect of noise. Edge gradient of these operators is lying in the 

range [0,2;r] . The second are compass operators, that measure gradient 
in pre-defined directions using templates, and assign the output and 
the associated direction of the mask giving maximal absolute value 
response as the gradient magnitude and orientation respectively. 
30 Larger masks enable to quantize the orientation to greater 
precision . 
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Alternative edge detection techniques such as zero crossing of 
second derivatives (such as a Laplace operator) may be considered 
also . 

5 

It should be addressed that the orientation computed by edge 
gradient measures refers to collinear orientation, i.e. if two 
orientations are collinear (that is, the difference between them is 
71, or equivalently , they have opposite edge sense), they should be 

10 considered as the same orientation. Therefore, in the absence of a 

model of the object (that e.g. would pair-wise group and label upper 
and lower rib edges in a thorax image) , which may be used to 
interpret the local orientation of a pixel and assign that pixel to 
its object boundary, the domain of the orientation 6 is [0 ? ^*]. For 

15 analyzing patterns with constant gradient orientation, this domain 
may be sufficient. For example, for analyzing the gray value 
trabecular texture pattern inside tubular bones such as the femur, 
which has a local preferential direction, or for analyzing the 
preferential direction of cortical (outer) edges of tubular bones, 

20 the edge gradient measure is useful for determining the main 
orientation . 

However, when determining one of the four orientations of a 
thorax (upside/down, left/right, and reversals), the edge 

25 orientation is insufficient when quantized to the domain [0,*] to 
discriminate between the upright and the upside-down direction or 
the left-to-right and the right-to-left orientations. 

Each rib is composed of a pair of (nearly parallel) edges; the 
orientation of juxtaposed edge pixels on either rib edge thus 

30 differs by 180 degrees. Hence, a histogram of quantized edge 
orientations will show two peaks, each one associated with a 
dominant edge orientation, and therefore it is impossible to 
discriminate between flipped or 180 degree rotated image versions. 

35 In the sequel, measures that have a domain of [0,2;r] are 

addressed, and that are insensitive to the sense of the edge. It 
means that, for example, any of a pair of juxtaposed thorax rib 
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edges will yield a similarly oriented analyzing shape feature 
vector, hence the orientation of the thorax can be discriminated 
into one of the four main orientations. 

The shape features outlined in the sequel are generalized to 3D 
space curves; the case of contours in the plane of a projection 
image (such as a radiography) is treated as a special case. 
Extensions to similar shape features of surfaces in 3D images (such 
as CT or MR) are also indicated. 

Curvature and torsion measure 



Arc length, curvature and torsion describe local properties of 
a general curve in 3D, the curvature of so-called principal normal 
sections describe the local properties of surfaces in 3D. The main 
15 tool for such investigations is a local coordinate frame, linked to 
a point on the curve or surface. 

Given an arbitrary parameterization 
y(t) = {x = (p x (t),y = (p 2 (t),z - (p 3 (t)} / and an associated natural 

20 parameterization y*(s) = {x = <p } (s),y = <p 2 (s),z = (p 3 (s)} of a space curve with 

s the arc length parameter, then the expressions for the tangent 

vector t , normal vector n and binormal vector b at a location p 0 

with natural parameter s 0 and arbitrary parameter t 0 are given by 



'(■So) 



n a - ~~~ t — t = b 0 x t 0 



Of\ — x — * II — n^. . . „ 



\\r'(t 0 )xr"(t 0 



Jointly t , n and b form a local coordinate system, called the 
30 Frenet-Serret frame. Tuples of these vectors define the osculating 

(Fand n) , normal (n and b) and rectifying plane (? and b). The 



orientation of the tangent t 0 is parallel to the tangent line to the 
curve through p 0 , and points in the direction defined by the sense 
in which the curve is traversed, a traversal in reverse sense 
yielding an opposite direction of the tangent vector. 

In the context of image processing, the direction of traversal 
of a curve may be based on the distribution of gray values across 
the edge, that is, for example such that the brighter image area is 
at the right hand side, and the darker image area is on the left 
side . 

The orientation of the normal n 0 (also called the main normal) 
is chosen such that the vector points in the direction where the 
projection of the curve on the osculating plane is concave. This 
orientation choice is independent of the traversal direction of the 
curve. The local Frenet -Serret frame at a point p Q of the curve is 
depicted in fig. 3 and the following relationships further hold on 
the depicted vectors : 

r'(t 0 )//f'(s 0 ) 

The change of the vectors t ,n and b when the point p 0 moves 
along the curve are given by the formulas of Frenet - Serret 
dl 

= KYI 

ds 

dn r 

— = -Kt + rb 
ds 

db 

— = -zn . 
ds 

The geometric meaning is that the parameter a: and r relate to 
the curvature and torsion of the curve respectively at point p 0 . 
They are given by the formulas, expressed in natural resp. 
arbitrary parameterization 
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H|f* ,, (*o)|| = 



\\r' (Qxf"(<o) 



T = 



(fXs 0 )xf'\s 0 ))-r\s 0 )_(fXto>f"(to))-r n Vo).. 



tt'Co) <p 2 '(t 0 ) <p 3 \t 0 ) 

nV 0 ) <p 2 v 0 ) nVo) 

<pr(t 0 ) <p 2 m <,t 0 ) <p 3 'v 0 ) 



K 



A 2 +B 2 +C z 



with 



A = 



<('o) <'('o) 
<P 2 \t 0 ) <p 2 '\t 0 ) 



,B = 



<P,\t 0 ) <p^(t 0 ) 



,C = 



Obviously, a: is a scalar chosen such that that the magnitude of 
the normal vector n is unity. 



The reciprocal of the curvature is p = — , the radius of 

K 

10 curvature. The curvature of a straight line is zero, while the 

curvature of a circle is constant and equal to the reciprocal of the 
circles' radius. 



From the formulas of Frenet-Serret , it can be seen that the 
15 curvature is equal to the angular velocity of the tangent vector, 
and hence is a measure of the rate of change of the tangent vector 
as the arc length varies. The torsion is equal to the angular 

velocity of the binormal vector b , and hence is a measure of the 
twisting of the curve out of the osculating plane. 

20 

Plane (2D) curvature computation 

For two-dimensional curves, in the plane of the image, the 
torsion is zero everywhere, and the curvature for a general 
parameterization reduces to 

|*W'(0-*'W(0| 
* (*'(0 2 +/(0 2 ) 3/2 ' 

(for a natural parameterization the denominator is 1) . In a 
digital image, derivatives of a spatially sampled version of the 
curve are estimated using finite differences. 
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Note that, although the triple (t,n,b) forms a right-handed 
coordinate system, and because /cis strictly positive, the plane 

coordinate system (t ,n) does not inherently imply a right-handed 
coordinate system. 

5 

As stated higher, the topology component of the application 
cannot be solved without knowing the characteristics of the shape of 
which the topology is requested. To determine the shape's topology, 
only pixels of the shape having curvature and torsion 

w characteristics within a specific spatial range may be taken into 
account. Moreover, the characteristics themselves, which must be 
fulfilled, must be chosen in relationship to the analyzing scale of 
the topology. On the one hand, topology of the shape at a too fine 
scale level may be considered irrelevant, because it is mostly 

15 attributable to noise (either intrinsic to the shape, or extrinsic 
such as image quantization noise) . On the other hand, topology of 
the shape at a too coarse level may also be unwanted, because it may 
lead to missing the discriminating topological characteristic. It is 
clear that when the suitable scale is selected, image areas having 

20 curvilinear structures that are more curved will contribute more to 
the determined image orientation. 



Fourier-based curvature estimation 
25 Curvature measures may also be computed from a complex contour 

representation. There are two basic approaches, depending on the 
specific representation 
o W (0 = (x(t),y(t)) 

If X(f)and Y(f) denote the Fourier transform of the contour 
30 coordinate signals x(0 and y(t) , using the Fourier derivative 

property that holds on a f ourier transform F(f) of a signal , 

namely F(f) = j27rfF(f),F(f) = -(27rf) 2 F(f), the plane equation of 
curvature can be rewritten as 

m = F-\JX(f))F-\f 2 Y(f))-F-\f 2 X(f))F-\fr(f)) 
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• k(0 = *(0 + /v(0 



Using the derivatives of u(t) , i.e. u(t) = x(f) + jy(t) and 
= Jc(7) + jy(£) , the plane equation of curvature can be rewritten as 



wherein the derivatives u(t) 9 u(t) are estimated from the inverse 
of the filtered Fourier transforms of the complex signal u(t) . 

Because contours are extracted from noisy images, and are 
spatially sampled, the high-pass filters that implement the 
derivatives accentuate the influence of high-frequency noise, which 
could undermine the curvature estimation completely. Applying multi- 
scale low-pass filters prior to estimating the derivatives can 
circumvent the noise-enhancing nature of the derivatives. Gaussian 
filtering, with the standard deviation parameter cr of the time- 
domain Gaussian is associated with the analyzing scale and controls 
the degree of contour smoothing (low-pass filtering) . However, the 
Gaussian filtering modifies the signal amplitudes of the filtered 
contour, implying the so-called shrinking effect. Hence it also 
affects estimated curvature, since the latter depends on the curve 
scale. The shrinking effect can be prevented by either re- 
normalizing the contour through the principle of spectral energy 
conservation, or through the principle of perimeter conservation. 

Preferential 2D image orientation based on normal vector and 
curvature computation 

Based on the foregoing differential geometry basics, 
preferential orientation of an image structure with respect to a 
given analyzing coordinate system (e.g. one that is aligned with the 
original image axes, i.e. Cartesian) may now be expressed by drawing 
the reference coordinate system locally, that is, attached to the 
point p 0 , and determining the quadrant of the analyzing coordinate 

system in which the normal n 0 lies. 



-lm{u(t)ii*(t)} 
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For example, assigning quadrants as given in fig. 6, a curve, 
which has everywhere convex shape with respect to the lower image 
border of the analyzing coordinate system, will have all normal 
vectors along its path either assigned to quadrant 3 or quadrant 4. 

5 Because quadrants 3 and 4 are associated with the negative y values, 
the overall orientation sense of the shape is determined to be in 
the direction of the lower image border. 

Likewise, all normal vectors of a convex shape with respect to 
the left image border, will have quadrant assignments 2 and 3. The 

w distribution of quadrant votes are subsequently collected in a 

histogram, and the quadrant or pair of quadrants with highest votes 
is taken to assign the principal shape orientation. To weight points 
with high curvature as more important indicators of the orientation 
sense of the shape, each vote may be weighted by the value of the 

15 curvature, which is an indicator of the orientation's strength. 

The direction of the center of the osculating circle, which is 
in the direction of the normal, may be determined by constructing a 
left and right chord to point p 0 , each with equal arc length 

20 distance to point p 0 , and performing the subtraction of these chord 
vectors (when the chord vectors are drawn both emanating from the 
point p 0 , the subtraction is performed as an addition) . The 
subtraction vector amounts to a second derivative (the chord length 
vectors being first derivatives) , and when the natural 

25 parameterization is used, the subtraction vector points in the 
direction of the normal vector. 

This vector operation is depicted in fig. 4(b). 

For a circle, which has constant curvature along the curve, and 
when a natural parameterisat ion is used, the chords associated with 
so points at equal arc length distance with respect to p Q , will have 
equal length. Therefore, their addition vector points in the 
direction of the circle's centre, which direction is perpendicular 
to the tangent vector. 

For a straight line, and when a natural parameterisat ion is 
35 used, addition of the chords yields the zero vector, which is 

plausible because a line has zero curvature everywhere. Either the 



GN03029 



2004-02-061:56 pm 



- 17 - 

length of the addition vector, or the angle subtended between the 
chords may be used as a curvature measure, with larger lengths or 
smaller angles respectively denoting higher curvature. The direction 
of the addition vector, corresponding to an associated subtended 

5 chords' angle of less than 180 degrees is used to determine the 

convex side of the curve, or equivalent ly is used to determine the 
normal direction, the sense of which is such that when looking 
towards the vector's end point, the curve is viewed concavely. This 
subtended angle 0 ik is depicted in fig. 4(a) and is computed on the 

10 basis of the ^-vectors (Jc referring to the arc length distance, 
expressed as a number of contour pixels, and determining the 
analyzing scale) 

a ik = (i f - - x i+t9 y i - y i+k ) 

K =(xt-xi-k>yi-yi-k) 
( - =* ^ 

is 9 ik = arccos 







I 5 * I' 





cos(9 ik ) will be 1 for the angle at a cusp (corresponding to 
infinite curvature) , 0 for a 90 -degree angle and -1 along a straight 
line . 



20 As stated before, in digital images, special attention has to 

be paid to the scale at which the derivatives in the curvature 
formula are evaluated. This scale may be controlled by the choice of 
the arc length differences defining the chord vectors. Care must be 
exercised in selecting the right analyzing scale. Obviously, when 

25 the arc length distance between the analyzing points is too small, 
curvature due to noisy variations of the curve is measured. 
Conversely, when the distance is too large, curvature-based shape 
features may still have poor discriminatory power. The 
differentiation acts as a high-pass filter that accentuates the 

30 influence of high frequency noise, which can undermine the curvature 
estimation. Applying larger inter-pixel arc length distances in the 
evaluation of the differences has a similar effect as applying a 
low-pass filter in the derivatives estimation step, and will yield 
more stable curvature estimates. 
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Preferential 3D space curve orientation 

Similarly to using curvature by assigning the normal vector to 
a local quadrant of a local analyzing 2D coordinate system for 
5 assessing the shape orientation of a curve in 2D, normal and 

binormal may be used to assign them to octants of a local analyzing 
3D coordinate system, for assessing the shape orientation of a 3D 

curve or surface. For example, at each point p 0 of a 3D space curve, 

a normal plane is constructed at p 0 , formed by the normal vector n Q 

w and binormal vector b 0 . The addition vector n 0 + b 0 lies in the normal 
plane and its orientation sense may be quantized into one of the 
eight analyzing Cartesian space quadrants. Analogous to the 2D 
case, the weighting may be generalized by multiplying the vote with 

a scalar depending on the local curvature and torsion at p Q . 

15 

Application of image orientation detection 
1. Thorax orientation detection 

Curvature of ribs or the ribcage may be used to infer the 
20 orientation and reflection of the thoracic X-ray image. The 

curvature is calculated for all edge segments resulting from an 
edge-based segmentation and the preferential direction (top-down, 
right-left, left-right or down-top) may be inferred on the basis of 
a voting scheme of analyzing coordinate system quadrant assignments. 
25 Edge segments associated with direct exposure contours will 

typically have less influence because their number is much smaller 
than edge segments belonging to ribs. Alternatively, they may be 
discarded based on their classification as being a diagnostic/direct 
exposure contour. 

30 

According to fig. 5, the algorithm for detection of the 
orientation in one of four classes (upright, upside-down, right- to- 
left with upper part of the thorax at the right side, left to right 
with upper part of the thorax at the left side) is as follows: 
35 Step 1 : segmentation of the lung fields using techniques of the 

prior art such as thresholding or region-growing (region-based 
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methods) , or active shape models (contour-based methods) . Only image 
data pertaining to the lung fields may be used in further process 
steps . 

Step 2 : edge detection using techniques of the prior art such 
5 as the Sobel operator, and thinning of edges by non-maximum 
suppression to yield one-pixel thick edges. 

Step 3 : edge -following to group touching edge pixels into edge 
segments . 

w Step 4 : curvature and normal computation at each edge point 

along the edge segment according to the principles and technique 
outlined before. 

Step 5 : quantizing the normal direction. This operation is 
equivalent to coordinate quadrant voting when there are 4 angular 

is intervals corresponding to partitioning of the image plane using 

axes parallel with the coordinate axes. In the example of fig. xxx, 
the normals pointing towards the center of the osculating circle, 
mostly point in the direction of quadrant 4, edge pixels nearer to 
the mediastinum have their normals pointing towards quadrant 3. 

20 Similar voting result holds true for the other ribs in this lung 
area; the reverse holds true for the associated ribs in the 
contralateral lung area. Quadrant voting is plausible when a 
topology modification is considered that will change the orientation 
of the image by a multiple of 90 degrees with respect to the axes of 

25 the current image coordinate system. 

Step 6 : grouping of angular interval entries to cover the 
desired image orientations. For thorax examinations, which normally 
are taken in the upright position (such that the thorax symmetry 
line is a vertical, and the upper thorax part is at the top of the 

30 image), the group sums of quadrant pairs 1+2, 2+3, 3+4, 4+1 are 

considered. Because quadrants 3 and 4 are below the x-axis, and are 
both touching to the lower image border (when the coordinate axes 
are drawn in the center of the image) , the group sum 3+4 is 
associated with the upright thorax position. The group votes 

35 according to quadrants 1+2 are associated with the upside-down 

thorax position; similarly, group votes of quadrants 2+3 and 4+1 are 
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associated with the right -to- left resp. left- to-right thorax 
position . 

Step 7 : selecting the image orientation according to the 
maximal group vote. According to fig. 5 the image is in the upright 
5 position when the group sum of quadrant 3 and 4 votes is maximal 
among the sum 1-2, 2-3, 3-4, 4-1. 

The orientation of a frontal and a lateral skull contour may 
also be inferred based on the cumulative direction of normal 

10 vectors. Because of the symmetrical shape of the skull contour, 

addition of all normal vectors belonging to the diagnostic/direct 
exposure contour, will have canceling components in the 
perpendicular direction, however, in the direction corresponding to 
the skull symmetry line, they will have a cumulative additive 

is effect. 

This principle may also be used to detect 3D preferential 
orientation of spherical shapes in 3D voxel data, such as occurring 
frequently in CT or MR cerebral examinations. 

20 Curvature and torsion are interesting descriptors for 3D space 

curves such as catheters, and the curve corresponding to the spinal 
center . 

Curvature and torsion may be used to characterize spinal 3D 
25 shape and to infer reflection of a lateral thorax view. A lateral 
view of the spinal chord has positive curvature for its thoracic 
part (the medical nomenclature of this curvature being kyphosis) , 
and has negative curvature for its lumbar part (the medical 
nomenclature of this curvature being lordosis) . Hence lateral views 
30 of the spine may be displayed or printed such that it has the 
standard or preferred orientation. 

From the foregoing, it will be clear that the disclosed 
algorithm to detect the orientation of thorax examinations is 
35 insensitive to the presence of collimation areas, for it uses edges 
as the main detected image feature. The influence of collimation 
borders is minimized for two basic reasons. First, the number of 
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edge pixels due to collimation borders is typically much less than 
the number of edge pixels due to edges. Hence their contribution in 
the voting is smaller. Secondly, the collimation borders are 
typically straight lines, meaning that their curvature is near to 
zero. Hence the weight attached to them is also very small. Hence, 
it is an advantage of the present invention that collimation need 
not be removed necessarily for correctly determining orientation. 
This feature is particularly advantageous in pediatric radiology, 
where the border may cover more than 3 0% of the image area. 



A second characteristic of the disclosed method is that it 
relies on image features which are very typical of a thorax image 
i.e. it relies on the shape of ribs, which are almost always present 
on a thorax, even when the contrast of some may partially be 
is diminished due to hyper-dense lung tissue. 

The algorithm further handles four types of thorax images: 
adult posterio-anterior (PA) and antero-posterior (AP) , and 
pediatric PA and AP projections. Although all of them are thorax 

20 images, their range of appearance differs significantly. 

Particularly pediatric images are most different from the standard 
adult thorax image. In the latter, the mediastinum is usually 
centered within the image, and the image itself contains only the 
thorax. In pediatric thorax acquisitions, the image may or may not 

25 include the head or part of it, the location of the arms is random, 
and part of the abdomen may sometimes be included and the area of 
exposed sub-diaphragm may differ. For these regions are relative 
free of the presence of structured edges, they will contribute 
negligibly to the voting measure. The algorithm is further robust 

30 against non-centered image and the presence of a centerline that may 
have a random angle with respect to the vertical axis of the image. 
Because the ribs are circular 3D structures spanning between the 
sternum and the spine, the curvature measure is also relative 
invariant to angular deviations of the centerline. 



35 



2. Mammographic orientation detection 
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Curvature of the mammographic skin line may be used to infer 
the orientation and reflection of the mammographic X-ray image. The 
curvature is calculated for all skin border edge segments resulting 
from an edge-based segmentation and, similarly to the thorax 
5 orientation detection, the preferential direction (top-down, right- 
left, left-right or down-top) may be inferred on the basis of a 
voting scheme of analyzing coordinate system quadrant assignments. 
The skin line edge can be determined by selecting and chaining the 
pixels on the region transition obtained by a segmentation such as 
w disclosed in European patent application 887 769. 

According to fig. 6, the algorithm for detection of the 
orientation of a mammographic CC (craniocaudal) view into one of 
four classes (corresponding to whether the chest wall is aligned 
is with the left, right, upper or lower image border) is as follows: 

Step 1 : segmentation of the skin line of the mammographic 
breast mass using techniques of the prior art such as thresholding 
or region-growing of the direct exposure area. 

Step 2 : edge detection to determine the one-pixel thick edges 
20 of the skin line (transition between breast mass and direct exposure 
area) . 

Step 3 : edge-following to group touching edge pixels into edge 
segments . 

Step 4 : curvature and normal computation at each edge point 
25 along the edge segment according to the principles and technique 
outlined before. 

Step 5 : quantizing the normal direction. This operation is 
equivalent to coordinate quadrant voting when there are 4 angular 
intervals corresponding to partitioning of the image plane using 
30 axes parallel with the coordinate axes. 

Step 6 : grouping of angular interval entries to cover the 
desired image orientations. 

Step 7 : selecting the image orientation according to the 
maximal group vote. In the example of fig. 6a, the normals around 
35 the breast nipple of the RCC view given mostly point in the 

direction of quadrant 1 and 4 (representing the right image border) . 
For the LCC view given, the normals around the breast nipple point 
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in the direction of quadrant 2 and 3 (representing the left image 
border) . For the RMLO and LMLO view given, these normals point 
towards quadrant 1 (representing the top-right corner) and quadrant 
2 (representing the top-left corner) respectively. 

5 

Variations of the outlined algorithm can be envisaged without 
departing from the scope of the present invention. For example, 
instead of computing the normal for all edge pixels and applying a 
voting mechanism, a single circle can be fitted to the edge data and 
w the direction of the midpoint of the fitted circle segment towards 
the circle center can be used as the main object orientation. This 
circle resembles a least squares approximating osculating circle of 
maximal curvature . 

15 Region-based shape orientation measures 

The former methods actually require that a curve be extracted 
from the image data, for example by edge detection. However, in the 
absence of a well defined boundary of a curved object, the curvature 
20 of the overall shape of the object can also be computed by using 
iso- intensity lines embedded in the object. Hence, the voting 
strategy outlined before, may now be applied to the group of pixels 
inside the object using the local quantized curvatures of the 
available or selected object pixels. 

25 

Region representation 

In its simplest form, a region may be viewed as a grouping or 
collection of pixels belonging to an entity having a problem- 
specific semantic (e.g. all pixels belonging to an object part). At 

30 a higher level of abstraction, a region may be described by its 
decomposition in smaller primitive forms (such as polygons, or 
quadtrees) . A region may also be described by its bounding region 
such as the Feret box, the minimum enclosing rectangle or the convex 
hull. Finally, a region may be represented by its internal features 

35 such as its skeleton, or a run-length representation. 

Iso- intensity curvature measure 
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An iso-intensity line is a curve of constant intensity, 
sometimes called an isophote. The geometric properties at a point 
(x,v)can be described by approximating the gray value function 
f(x,y)by its Taylor series 

Ar n2 

f(x + Ax, y + Ay) = f(x, y) + 
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Defining the gradient of the gray value function f(x,y) as 

m*,y) = 



~f( x >y) 




'fx' 






fy. 



the first derivative can be computed in any direction v as 
/ V =v-V/. 

10 The partial derivatives in the foregoing and following 

expressions are computed by convolving the original intensity image 
with directional derivatives of a Gaussian at a certain scale 
represented by the standard deviation cr of the Gaussian, that is 
based on the equality 

is 4- (/(*• y) ® G <*> y> a )) s y"> ® (iL G ( x > y> a ^ 

dx y {dx J 

The spatial extent of the Gaussian has a smoothing effect, and 
is advantageous to ignore curved structures that are attributable to 
noise or irrelevant image detail such as local trabecular bone 
structure . 

20 The normal to the isophote in a pixel is the gradient 

direction, which direction is tangent to the flowline through that 
pixel. The tangent to the isophote curve is the gradient 
perpendicular. The unity normalized vectors corresponding to these 

directions are given by w = n = Vf /|[ V/||, v =t=n ± . The orientation or 

25 sense of the normal n is chosen such that the vector points in the 
direction where the projection of the isophote curve on the tangent 
line is concave. This orientation choice is independent of the 
traversal direction of the curve. 

By taking the second-order term in the local Taylor 

30 approximation at a point (x, y) into account, second-order geometric 
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20 



H f (x,y) = 



properties can be computed of the local image gray value surface. 
The second-order term is described by 
the symmetric Hessian H f (x,y) 

fixity) f xy ( x >yY 

f x y(X*y) fyy(x>y)_ 

This matrix can be used to compute the second derivatives in 
other directions, such as the local gradient w , tangent v 
coordinate axes, as follows 

f m = v T Hw , f n = v T Hv , f m = w T Hw 

These second-order derivatives are the constituent elements of 
the Hessian matrix in the vw coordinate system, obtained by applying 
the rotation matrix R to the xy coordinates, 
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RH f R T = 



1 



' f f 

f f 

" flfyy " 2f x f y f xy + (/, 2 - //) + fJ,(J m - f„) 

(fy -/?) + /,/,(/« - /„) /x 2 /„ + *fJyf* + flfyy _ 
The Hessian is real and symmetric and has the properties that 
its determinant is equal to the product of its eigenvalues, and is 
invariant to the selection of the original axes xandy . 

For v points in the isophote contour direction, the second 
derivative in the contour direction SDCD is given by 

f 2 f 

J x J yy 

vv 



fx f vv ^fxfyfxy fy fxx 
fl+fl 



The second derivative in the isophote contour direction 
v = (f y ,—f x ) is related to the isophote curvature defined as the change 
in orientation of the tangent angle as one moves along the isophote 
path 
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where yjjf + f* is the gradient magnitude. 



10 



Similar to the voting procedure outlined in the curve-based 
orientation detection method, the orientation direction of a curved 
grey value structure is computed by a weighted voting of quantized 
local orientations of the directions corresponding to the normal 
vector. The weight attached to the vote is proportional to the value 
of the isophote curvature K , meaning that image areas having iso- 
intensity lines that are more curved will contribute more to the 
determined image orientation. 



15 



20 



25 



This method is extended to determine the main orientation of a 
3D iso-grey surface patch in a 3D image. At each point (x,y,z) the 
principal curvatures K[,?q {fq>iq) and their corresponding surface 
vectors c, and c 2 lying in the plane tangent to the iso-grey surface 
are computed as follows. First, the gradient vector g = (f x ,f y9 f z ) is 
calculated, which vector is aligned with the surface normal vector, 
and the Hessian matrix 



H = 



f f f 

J xx J xy J xz 

f f f 

J xy J yy J yz 

f xz fyz f zz 



with derivatives evaluated by convolving the image with 
derivates of Gaussian kernels. Next, the Hessian is rotated to align 
the first axis with the gradient. The resulting matrix has the form 



H' = 



0 f m f u 
0 /„ / W J 



0 



ss 

0 H T 



with /the second derivative in the gradient direction and 
H T a 2D Hessian in the touching plane perpendicular to the gradient. 
Finally, the eigenvalues and eigenvectors of the Hessian H T are used 
to describe second-order geometry of the local gray value surface of 
the 3D image. The first eigenvector ^corresponding to the largest 
absolute eigenvalue X is the direction of the greatest curvature 
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(second derivative). The second eigenvector^, whose corresponding 
eigenvalue has the smallest absolute value, is the direction of 
least curvature. These eigenvectors lie in the touching plane 
perpendicular to the gradient. The corresponding eigenvalues are the 

5 respective amounts of these curvatures, denoted by ^respectively, 
and referred to as principal curvatures. They are invariant under 
rotation of the original coordinate system. The eigenvectors are 
called principal directions p ? gand are perpendicular to each other. 
They are directions of pure curvature for their mixed partial 

10 derivatives are zero. For curvature is defined along a line, the 
curvatures /q,/^ of the intersection curves C,andC 2 of the surface 
patch with the planes (c l9 g) and ic 2 ,g) respectively are related as 
in 2D with the second derivatives in contour direction. 
_ \ SDC X D 

K ' = ~Uf~ III 

/L SDC.D 



is The quantity represented by the determinant of the Hessian// r 

is called the Gaussian curvature. The trace of the Hessian (which is 
the sum of the diagonal elements) is also invariant to the selection 
of xandy . Half the trace is equal to the mean of the principal 

curvatures /q,/^ , and are invariant to the selection of xand^ • 
20 Similar to the voting procedure outlined in the curve-based 

orientation detection method, the orientation direction of a curved 
3D grey value structure is computed by a principal curvature- 
weighted voting of quantized local orientations vectors. For 
curvature is defined along a line, the curvature value used at the 
25 vote weight is a combination of the principal curvatures. Suitable 
combinations, rewritten in term of derivatives with respect to the 

principal directions /?and# are the Gaussian curvature K = k x k 2 = L pp L qq 

(which geometrically represents the extra amount of area enclosed by 
the perimeter of the curved surface patch w.r.t. the area of a flat 

30 patch), the mean curvature H = (L +L )/2 , or the Laplacian 
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L = L +L ~2H . They are invariant under rotation of the coordinate 
system. 
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